
//extraterrestrial radiant flux
function get_e0(day_number){

    var e0;

    e0 = 1367*(1+0.033*Math.cos(360*(day_number-3)*d_to_r/365));

    return e0;
}



 //part of year
 function get_g(day_number){
    var g;
    g = 360*((day_number-1)/365);
    return g;
}

//equation of time
function get_eT(g){
    var eT = 2.2918*(0.0075+0.1868*Math.cos(g*d_to_r)-3.2077*Math.sin(g*d_to_r)-1.4615*Math.cos(2*g*d_to_r)-4.089*Math.sin(2*g*d_to_r));
    return eT;
}
//apparent solar time
function get_aST(lST, longitude, eT,zone) {
    var aST;
    if (zone*1<0)
{ aST=lST*1+eT/60+(-zone*15+longitude*1)/15;
} 
else
{ aST=lST*1*1+eT/60+(longitude-zone*15)/15;
}
return aST;
}

//declination
function get_delta(day_number) {
    var delta = 23.45*Math.sin((360*(day_number+284)/365)*d_to_r);
    return delta;
}

// Hour angle
function get_h(aST) {
    var h = 15*(aST-12);
    return h;
}

// solar altitude
function get_beta(latitude, delta, h) {
    var beta = (Math.asin(Math.cos(latitude*d_to_r)*Math.cos(delta*d_to_r)*Math.cos(h*d_to_r)+Math.sin(latitude*d_to_r)*Math.sin(delta*d_to_r)))*r_to_d;
    return beta;
}

// solar azimuth
function get_fi(h, delta, beta,latitude,aST) {
    var fi;
    var fi2;
    var fi4;
    var fi1 = (Math.asin((Math.sin(h*d_to_r))*(Math.cos(delta*d_to_r))/(Math.cos(beta*d_to_r))))*r_to_d;
    if (fi1<0) {fi2 = fi1 + 180;}
    if (fi1>0) {fi2 = fi1 - 180;}
    var fi3 = (Math.acos((Math.cos(h*d_to_r) * Math.cos(delta*d_to_r) * Math.sin(latitude*d_to_r) - Math.sin(delta*d_to_r) * Math.cos(latitude*d_to_r))/Math.cos(beta*d_to_r)))*r_to_d;
    fi4 = -fi3;
    if (fi1.toFixed(2) == fi3.toFixed(2)) {fi = fi1;} else {
    if (fi1.toFixed(2) == fi4.toFixed(2)) {fi = fi1;} else {
    if (fi2.toFixed(2) == fi3.toFixed(2)) {fi = fi2;} else {
    if (fi2.toFixed(2) == fi4.toFixed(2)) {fi = fi2;}
    }}}
if (aST<12 & fi>0) {fi=-fi;}
if (aST>12 & fi<0) {fi=-fi;}
    return fi;
}

/*Surface-solar azimuth. The surface azimuth ? is
defined as the displacement from south of the projection, on the horizontal
plane, of the normal to the surface. Surfaces that face west
have a positive surface azimuth; those that face east have a negative
surface azimuth */
function get_gamma(fi,psi) {
    var gamma = fi-psi;
    return gamma;
}

//Incidence angle
function get_Fi(beta,gamma,tilt) {
    var Fi = (Math.acos(Math.cos(beta*d_to_r)*Math.cos(gamma*d_to_r)*Math.sin(tilt*d_to_r)+Math.sin(beta*d_to_r)*Math.cos(tilt*d_to_r)))*r_to_d;
    return Fi;
}

// relative air mass
function get_m(beta){
    var m;
    if (beta>0) {
    m = 1 / (Math.sin(beta*d_to_r)+0.50572*Math.pow((6.07995+beta), -1.6364));
    }
    else {
        m = Infinity;//may be unlimited
    }
    return m;
}

//beam air mass exponent
function get_ab(taub,taud){
    var ab = 1.219 - 0.043*taub-0.151*taud-0.204*taub*taud;
    return ab;
}

//beam normal irradince
function get_eb(e0,taub,m,ab){
    var eb = e0*Math.exp(-taub*Math.pow(m, ab));
    return eb;
}

//Surface beam irradiance
function get_etb(eb,Fi){
    var etb;
    if (Fi > 90 || Fi < -90 ) 
    { 
    etb = 0;
    }
    else {
        etb = eb * Math.cos(Fi*d_to_r);
    }
    return etb;
}
        
//Ratio Y of sky diffuse radiation on vertical surface to sky diffuse radiation on gorizontal surface
function get_y(Fi){
    var Y = Math.max((0.55 + 0.437 * Math.cos(Fi*d_to_r) + 0.313 * Math.cos(Fi*d_to_r)*Math.cos(Fi*d_to_r)), 0.45);
    return Y;
}

//Difuse air mass exponent
function get_ad(taub,taud){
    var ad = 0.202 + 0.852 * taub - 0.007 * taud - 0.357 *taud * taub;
    return ad;
}

//diffuse irradiance on gorizontal surface
function get_ed(e0,taud,m,ad,beta){
    var ed;
    if (beta<0) 
    { 
    ed = 0;
    }
    else
    {
        ed = e0*Math.exp(-taud*Math.pow(m, ad));
    }
    return ed;
}

//diffuse irradiance on vertical surface
function get_etd(ed,y,tilt){
    var etd;
    if (tilt <= 90) 
    { 
    etd = ed * (y * Math.sin(tilt*d_to_r) + Math.cos(tilt*d_to_r)); 
    }
    else
    { 
    etd = ed * y * Math.sin(tilt * d_to_r);
    }
    return etd;
}
//ground reflected irradiance
function get_etr(eb,beta,ed,ground_reflect,tilt){
    var etr = (eb * Math.sin(beta*d_to_r) + ed) *  ground_reflect * (1- Math.cos(tilt * d_to_r))/2;
    return etr;
}

//Subtotal diffuse
function get_estd(etd,etr) {
var estd = etd + etr;
return estd;
}
    
//total surface irradiance
function get_et(etb,etd,etr){
    var et = etb + etd + etr;
    return et;
}

//sol-air temperature
function get_te(t0,alfa_h0,et,epsilon_deltaR_h0){
    var te = t0*1 + alfa_h0 * et - epsilon_deltaR_h0;
    return te;
}

//Calculate sun position and incidence angle
function sun_on_plane(day_number, lST, latitude, longitude,zone, psi,tilt) {  //lST - local standart hour

this.e0 = get_e0(day_number);  //extraterrestrial radiant flux
this.g = get_g(day_number); //part of year
this.eT= get_eT(this.g); //equation of time
this.aST = get_aST(lST, longitude, this.eT, zone); //apparent solar time
this.delta = get_delta(day_number); //declination
this.h = get_h(this.aST); // Hour angle
this.beta = get_beta(latitude, this.delta, this.h); // solar altitude
this.fi = get_fi(this.h,this.delta,this.beta,latitude,this.aST); // solar azimuth
this.gamma = get_gamma(this.fi, psi); //Surface-solar azimuth
this.Fi = get_Fi(this.beta,this.gamma,tilt); //Incidence angle
}


//Calculate surface irradiance and sol-air temperature
function radiation(day_number, lST, latitude, longitude,zone, psi,tilt,taub,taud,ground_reflect,t0,alfa_h0,epsilon_deltaR_h0) {
    mysun_on_plane = new sun_on_plane(day_number, lST, latitude, longitude,zone, psi,tilt);//Create sun_on_plane
    this.m = get_m(mysun_on_plane.beta); // relative air mass
    this.ab = get_ab(taub,taud); //beam air mass exponent
    this.eb = get_eb(mysun_on_plane.e0,taub,this.m,this.ab);//beam normal irradince
    this.etb = get_etb(this.eb,mysun_on_plane.Fi);//Surface beam irradiance
    this.y = get_y(mysun_on_plane.Fi);//Ratio Y of sky diffuseradiation on vertical surface to sky diffuse radiation on gorizontal surface
    this.ad = get_ad(taub,taud);//Difuse air mass exponent
    this.ed = get_ed(mysun_on_plane.e0,taud,this.m,this.ad,mysun_on_plane.beta);//diffuse irradiance on gorizontal surface
    this.etd = get_etd(this.ed,this.y,tilt);//diffuse irradiance on vertical surface
    this.etr = get_etr(this.eb,mysun_on_plane.beta,this.ed,ground_reflect,tilt);//ground reflected irradiance
	this.estd = get_estd(this.etd,this.etr);//Subtotal diffuse
    this.et = get_et(this.etb,this.etd,this.etr);//total surface irradiance
    this.te = get_te(t0,alfa_h0,this.et,epsilon_deltaR_h0);//sol-air temperature
}

